module tracer_mod

	use const_mod
	use namelist_mod
	use mesh_mod
	use allocator_mod
	use parallel_mod


	implicit none

	private
	public tracer_type
	public tracer

	type tracer_type
		type(mesh_type), pointer :: mesh => null()
		real(r8), allocatable, dimension(:,:,:)   :: u
		real(r8), allocatable, dimension(:,:,:)   :: v
		real(r8), allocatable, dimension(:,:,:)   :: we
		real(r8), allocatable, dimension(:,:,:)   :: we_exp ! explicit velocity
		real(r8), allocatable, dimension(:,:,:)   :: we_imp ! implicit velocity
		real(r8), allocatable, dimension(:,:,:,:) :: q
		real(r8), allocatable, dimension(:,:)     :: qanl       ! exact solution, only for 2D stationary vertex and 2d moving vertex test
		real(r8), allocatable, dimension(:,:,:)   :: qmf_lon
		real(r8), allocatable, dimension(:,:,:)   :: qmf_lat
		real(r8), allocatable, dimension(:,:,:)   :: qmf_lev
		real(r8), allocatable, dimension(:,:,:)   :: divx
		real(r8), allocatable, dimension(:,:,:)   :: divy
		real(r8), allocatable, dimension(:,:,:)   :: cflx
		real(r8), allocatable, dimension(:,:,:)   :: cfly
		real(r8), allocatable, dimension(:,:,:)   :: cflz
		real(r8), allocatable, dimension(:,:,:)   :: qx
		real(r8), allocatable, dimension(:,:,:)   :: qy
		real(r8), allocatable, dimension(:,:,:)   :: mfx
		real(r8), allocatable, dimension(:,:,:)   :: mfy
		real(r8), allocatable, dimension(:    )   :: tm ! total mass
		real(r8), allocatable, dimension(:,:,:)   :: fx_low ! horizontal FCT
		real(r8), allocatable, dimension(:,:,:)   :: fy_low
		real(r8), allocatable, dimension(:,:,:)   :: fz_low ! vertical FCT 
		real(r8), allocatable, dimension(:,:,:)   :: q_td ! horizontal and vertical FCT
		real(r8), allocatable, dimension(:,:,:)   :: rin
		real(r8), allocatable, dimension(:,:,:)   :: rout
  contains
  	procedure :: init  => tracer_init
  	procedure :: clear => tracer_clear
  	final :: tracer_final
	end type tracer_type
  
  type(tracer_type), allocatable :: tracer(:)

contains
  subroutine tracer_init(this, mesh)

		class(tracer_type), intent(inout) :: this
		type(mesh_type), intent(in) :: mesh
		integer itracer

		call this%clear()

		call allocate_array(mesh, this%u      , half_lon=.true., full_lat=.true., full_lev=.true.)
		call allocate_array(mesh, this%v      , full_lon=.true., half_lat=.true., full_lev=.true.)
		call allocate_array(mesh, this%we     , full_lon=.true., full_lat=.true., half_lev=.true.)
		call allocate_array(mesh, this%we_exp , full_lon=.true., full_lat=.true., half_lev=.true.)
		call allocate_array(mesh, this%we_imp , full_lon=.true., full_lat=.true., half_lev=.true.)

		allocate(this%q(mesh%full_lon_lb:mesh%full_lon_ub, mesh%full_lat_lb:mesh%full_lat_ub, mesh%full_lev_lb:mesh%full_lev_ub, ntracers))
		allocate(this%qanl(mesh%full_lon_lb:mesh%full_lon_ub, mesh%full_lat_lb:mesh%full_lat_ub))
		call allocate_array(mesh, this%qmf_lon, half_lon=.true., full_lat=.true., full_lev=.true.)
		call allocate_array(mesh, this%qmf_lat, full_lon=.true., half_lat=.true., full_lev=.true.)
	  call allocate_array(mesh, this%qmf_lev, full_lon=.true., full_lat=.true., half_lev=.true.)
	  call allocate_array(mesh, this%divx   , full_lon=.true., full_lat=.true., full_lev=.true.)
	  call allocate_array(mesh, this%divy   , full_lon=.true., full_lat=.true., full_lev=.true.)
	  call allocate_array(mesh, this%cflx   , half_lon=.true., full_lat=.true., full_lev=.true.)
	  call allocate_array(mesh, this%cfly   , full_lon=.true., half_lat=.true., full_lev=.true.)
	  call allocate_array(mesh, this%cflz   , full_lon=.true., full_lat=.true., half_lev=.true.)
	  call allocate_array(mesh, this%qx     , full_lon=.true., full_lat=.true., full_lev=.true.)
	  call allocate_array(mesh, this%qy     , full_lon=.true., full_lat=.true., full_lev=.true.)
	  call allocate_array(mesh, this%mfx    , half_lon=.true., full_lat=.true., full_lev=.true.)
	  call allocate_array(mesh, this%mfy    , full_lon=.true., half_lat=.true., full_lev=.true.)

	  allocate(this%tm(ntracers))

	  if (upwind_hori_fct) then
	  	call allocate_array(mesh, this%fx_low , half_lon=.true., full_lat=.true., full_lev=.true.)
	  	call allocate_array(mesh, this%fy_low , full_lon=.true., half_lat=.true., full_lev=.true.)
	  end if
	  if (upwind_vert_fct) then
	  	call allocate_array(mesh, this%fz_low , full_lon=.true., full_lat=.true., half_lev=.true.)
	  end if
	  if (upwind_hori_fct .or. upwind_vert_fct) then
	  	call allocate_array(mesh, this%q_td   , full_lon=.true., full_lat=.true., full_lev=.true.)
	  	call allocate_array(mesh, this%rin    , full_lon=.true., full_lat=.true., full_lev=.true.)
	  	call allocate_array(mesh, this%rout   , full_lon=.true., full_lat=.true., full_lev=.true.)
	  end if

	end subroutine tracer_init

	subroutine tracer_clear(this)

		class(tracer_type), intent(inout) :: this

		if(allocated(this%u       )) deallocate(this%u      )
		if(allocated(this%v       )) deallocate(this%v      )
		if(allocated(this%we      )) deallocate(this%we     )
		if(allocated(this%we_exp  )) deallocate(this%we_exp )
		if(allocated(this%we_imp  )) deallocate(this%we_imp )
		if(allocated(this%q       )) deallocate(this%q      )
		if(allocated(this%qanl    )) deallocate(this%qanl   )
		if(allocated(this%qmf_lon )) deallocate(this%qmf_lon)
		if(allocated(this%qmf_lat )) deallocate(this%qmf_lat)
		if(allocated(this%qmf_lev )) deallocate(this%qmf_lev)
		if(allocated(this%divx    )) deallocate(this%divx   )
		if(allocated(this%divy    )) deallocate(this%divy   )
		if(allocated(this%cflx    )) deallocate(this%cflx   )
		if(allocated(this%cfly    )) deallocate(this%cfly   )
		if(allocated(this%cflz    )) deallocate(this%cflz   )
		if(allocated(this%qx      )) deallocate(this%qx     )
		if(allocated(this%qy      )) deallocate(this%qy     )
		if(allocated(this%mfx     )) deallocate(this%mfx    )
		if(allocated(this%mfy     )) deallocate(this%mfy    )

		if(allocated(this%tm      )) deallocate(this%tm     )

		if(allocated(this%fx_low  )) deallocate(this%fx_low )
		if(allocated(this%fy_low  )) deallocate(this%fy_low )
		if(allocated(this%fz_low  )) deallocate(this%fz_low )
		if(allocated(this%q_td    )) deallocate(this%q_td   )
		if(allocated(this%rin     )) deallocate(this%rin    )
		if(allocated(this%rout    )) deallocate(this%rout   )

	end subroutine tracer_clear

	subroutine tracer_final(this)
	  
	  type(tracer_type), intent(inout) :: this

	  call this%clear()
	    
	end subroutine tracer_final

end module tracer_mod